All external data files used in this code are available at https://github.com/MoTrPAC/motrpac-rat-training-mitochondria/tree/main/data
This script assumes that these files were downloaded to the data_dir path specified above
# read the mtDNA data file
mt_data = read.csv(
file=paste0(data_dir,"mtDNA_motrpac_pass1b.csv"),
stringsAsFactors = F
)
mt_data$bid = mt_data$Label
# read the RNA-seq metrics to obtain mtRNA read percentages
rnaseq_qa_qc = TRNSCRPT_META
pct_mito_reads = rnaseq_qa_qc$pct_chrM
names(pct_mito_reads) = as.character(rnaseq_qa_qc$vial_label)
rownames(rnaseq_qa_qc) = as.character(rnaseq_qa_qc$vial_label)
colnames(rnaseq_qa_qc) = tolower(colnames(rnaseq_qa_qc))
# load the pheno data
pheno = PHENO
cardiolipins_dir = paste0(data_dir,"/cardiolipins/")
cl_files = list.files(cardiolipins_dir,full.names = T)
cl_files = cl_files[!grepl("name_mapping",cl_files)]
cl_data = c()
cl_fdata = c()
for(f in cl_files){
f_arr = strsplit(f,split="_")[[1]]
tissue = f_arr[length(f_arr)]
tissue = gsub(".csv","",tissue)
tissue = toupper(tissue)
if(tissue=="MUSCLE"){
tissue = "SKMGN"
}
if(tissue=="hippocampus"){
tissue = "HIPPOC"
}
d = fread(f,stringsAsFactors = F,data.table = F,header=T)
d = d[!is.na(d[,1]),]
if(ncol(d)<51){next} # ignore inomplete tissues
# use bid instead of vial labels
colnames(d)[-1] = substring(colnames(d)[-1],1,5)
feature_d = cbind(tissue,d[,1])
if(length(cl_data)==0){
cl_data = d
cl_fdata = feature_d
}else{
cl_data = rbind(cl_data,d[,colnames(cl_data)])
cl_fdata = rbind(cl_fdata,feature_d)
}
}
cl_data = cl_data[,-1]
colnames(cl_fdata) = c("tissue","feature_ID")
cl_fdata = data.frame(cl_fdata,stringsAsFactors = F)
# transform WATSC and SKMGN to WAT-SC and SKM-GN
# to be consistent with the pheno data and the other
# data objects:
cl_fdata$tissue[cl_fdata$tissue=="SKMGN"] = "SKM-GN"
cl_fdata$tissue[cl_fdata$tissue=="WATSC"] = "WAT-SC"
cl_fdata$tissue[cl_fdata$tissue=="WAT"] = "WAT-SC"
table(cl_fdata$tissue)
##
## HEART KIDNEY LIVER LUNG SKM-GN WAT-SC
## 11 16 9 3 11 4
sort(table(cl_fdata$feature_ID),decreasing=T)[1:4]
##
## CL(72:8)>CL(18:2/18:2/18:2/18:2) CL(74:11)
## 4 4
## CL(74:9) CL(70:7)
## 4 3
cl_fdata[cl_fdata$feature_ID=="CL(72:8)",]
## tissue feature_ID
## 43 SKM-GN CL(72:8)
## 53 WAT-SC CL(72:8)
selected_cls = c(
"CL(72:8)>CL(18:2/18:2/18:2/18:2)",
"CL(72:8)",
"CL(74:11)",
"CL(74:9)"
)
# # Optional: get cardiolipins from the main landscape paper:
# sample_data = METAB_NORM_DATA_FLAT
# sample_data_featureinfo = sample_data[,1:5]
# # Limit the data to cardiolipins
# cl_inds = grepl("^CL",sample_data$feature_ID,perl = T)
# cl_data = sample_data[cl_inds,]
# cl_fdata = sample_data_featureinfo[cl_inds,]
# # change column names to bids
# m = unique(pheno[,c("pid","bid")])
# p2b = as.character(m[,2])
# names(p2b) = as.character(m[,1])
# colnames(cl_data)[-c(1:5)] = unname(p2b[colnames(cl_data)[-c(1:5)]])
# Add mito gene annotations
mt = fread(paste0(data_dir,
"external-datasets_gene-sets_Mitocarta_Rat.MitoCarta3.0.genes_only.txt"),
stringsAsFactors = F,data.table = F)
# additional pathway enrichment with mitocarta pathways
mt_pw = fread(paste0(data_dir,"Human.MitoCarta3.0.txt"),header = T,
stringsAsFactors = F,data.table = F)
# map transcripts to genes using biomart
library(biomaRt)
mart_rat = useMart("ensembl", dataset = "rnorvegicus_gene_ensembl")
attrs = c("ensembl_gene_id","ensembl_peptide_id","rgd_symbol")
rat_prot_gene = getBM(attributes = attrs, mart = mart_rat, uniqueRows=T)
rat_prot_gene = rat_prot_gene[rat_prot_gene$ensembl_peptide_id!="",]
rat_prot_gene = rat_prot_gene[rat_prot_gene$rgd_symbol != "",]
# read in rat to human mapping
rat_to_human = fread(paste0(data_dir,"RGD_ORTHOLOGS_20201001.txt"),
stringsAsFactors = F,data.table = F)
rat_to_human = rat_to_human[!is.na(rat_to_human$HUMAN_ORTHOLOG_SYMBOL),]
rat_to_human_map = unique(rat_to_human[,c("RAT_GENE_SYMBOL", "HUMAN_ORTHOLOG_SYMBOL")])
# make a pathway:member list
mitocarta_pathways = list()
for(pw in mt_pw$MitoPathway){
members = unname(unlist(strsplit(mt_pw[mt_pw$MitoPathway==pw, "Genes"], ', ')))
rat_members = rat_to_human_map[rat_to_human_map$HUMAN_ORTHOLOG_SYMBOL %in% members, 1]
rat_members = unique(na.omit(rat_members))
mitocarta_pathways[[pw]] = rat_members
#print(paste(pw, length(members), length(rat_members)))
}
save(rat_prot_gene,
mt,mt_pw,mitocarta_pathways,rat_to_human,
file=paste0(data_dir,"mito_gene_annotation.RData"))
Load above data:
load(paste0(data_dir,"mito_gene_annotation.RData"))
control_vials = rownames(pheno)[pheno$intervention=="control"]
control_vials = intersect(control_vials,names(pct_mito_reads))
baseline_pct_mito = data.frame(
pct_mito_reads = pct_mito_reads[control_vials],
tissue = pheno[control_vials,"tissue"],
sex = pheno[control_vials,"sex"]
)
baseline_pct_mito_male = baseline_pct_mito[baseline_pct_mito$sex == "male",]
baseline_pct_mito_female = baseline_pct_mito[baseline_pct_mito$sex == "female",]
currdf = baseline_pct_mito_male
tissue_med = tapply(currdf$pct_mito_reads,currdf$tissue,median)
tissue_med = sort(tissue_med,decreasing = T)
currdf$color = TISSUE_COLORS[currdf$tissue]
currdf$tissue = factor(currdf$tissue,levels = names(tissue_med))
currdf_male = currdf
currdf_male$sex = "Male"
currdf = baseline_pct_mito_female
tissue_med = tapply(currdf$pct_mito_reads,currdf$tissue,median)
tissue_med = sort(tissue_med,decreasing = T)
currdf$color = TISSUE_COLORS[currdf$tissue]
currdf$tissue = factor(currdf$tissue,levels = names(tissue_med))
currdf_female = currdf
currdf_female$sex = "Female"
baseline_pct_mito = rbind(currdf_male,currdf_female)
ggplot(baseline_pct_mito, aes(x=tissue, y=pct_mito_reads,fill=as.character(tissue))) +
geom_boxplot() + facet_wrap(~sex) + xlab("") +
theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) +
scale_fill_manual(
labels = sort(as.character(unique(baseline_pct_mito$tissue))),
values = TISSUE_COLORS[sort(as.character(unique(baseline_pct_mito$tissue)))])
mt_dna_copy = mt_data[mt_data$Group=="Ctrl",]
mt_dna_copy$sex = ifelse(mt_dna_copy$sex == "M","Male","Femlae")
tissue_med = tapply(mt_dna_copy$Level,mt_dna_copy$Tissue,median)
tissue_med = sort(tissue_med,decreasing = T)
mt_dna_copy$Tissue = factor(mt_dna_copy$Tissue,levels = names(tissue_med))
ggplot(mt_dna_copy, aes(x=Tissue, y=Level,fill=as.character(Tissue))) +
geom_boxplot() + facet_wrap(~sex) + xlab("") +
theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) +
scale_fill_manual(
labels = sort(as.character(unique(mt_dna_copy$Tissue))),
values = TISSUE_COLORS[sort(as.character(unique(mt_dna_copy$Tissue)))])
# Get data frames with all mito biomarkers, by tissue
mt_biomarkers = list()
for(tissue_name in unique(pheno$tissue)){
if(is.na(tissue_name)){next}
if(grepl("VENA",tissue_name)){next}
print(paste("Analyzing data from:",tissue_name))
# get all vial labels
tissue_vials = rownames(pheno)[pheno$tissue == tissue_name]
# limit to vial labels that appear in the rna-seq data
tissue_vials = intersect(names(pct_mito_reads),tissue_vials)
tissue_vials = setdiff(tissue_vials,OUTLIERS$viallabel)
if(length(tissue_vials)<5){next}
# get the bids
tissue_bids = sapply(tissue_vials,substring,1,5)
# save results in a data frame
tissue_df = data.frame(
mt_reads = pct_mito_reads[tissue_vials],
viallabel = tissue_vials,
bid = tissue_bids
)
tissue_df$mtdna = NA
if(sum(mt_data$Tissue==tissue_name)>0){
currmt_data = mt_data[mt_data$Tissue==tissue_name,]
rownames(currmt_data) = as.character(currmt_data$bid)
curr_shared_bids = intersect(tissue_df$bid,rownames(currmt_data))
for(bid in curr_shared_bids){
tissue_df[tissue_df$bid==bid,"mtdna"] = currmt_data[bid,"Level"]
}
tissue_df$mtdna = log2(tissue_df$mtdna)
}
tissue_df$clPC1 = NA
tissue_df$clPC2 = NA
tissue_df$clPC3 = NA
for(cl in selected_cls){tissue_df[[cl]] = NA}
curr_cl_data = NULL
if(sum(cl_fdata$tissue == tissue_name)>0){
curr_shared_bids = intersect(tissue_df$bid,colnames(cl_data))
curr_cl_data = cl_data[cl_fdata$tissue == tissue_name,curr_shared_bids]
metab_names = cl_fdata$feature_ID[cl_fdata$tissue == tissue_name]
curr_cl_data = aggregate(curr_cl_data,by=list(metab_names),mean,na.rm=T)
rownames(curr_cl_data) = curr_cl_data[,1]
curr_cl_data = curr_cl_data[,-1]
# remove columns that are all NAs
bid_na_rate = colSums(is.na(curr_cl_data)) / nrow(curr_cl_data)
curr_cl_data = curr_cl_data[,bid_na_rate<1]
curr_shared_bids = colnames(curr_cl_data)
# check missing values
na_rate = sum(is.na(curr_cl_data))/(nrow(curr_cl_data)*ncol(curr_cl_data))
if(na_rate >= 0.05){
curr_cl_data = NULL
}
else{
# use min imputation
curr_cl_data = t(curr_cl_data)
curr_cl_data[is.nan(curr_cl_data)] = NA
curr_cl_data[is.na(curr_cl_data)] = min(curr_cl_data,na.rm=T)
}
if(!is.null(curr_cl_data) && nrow(curr_cl_data)>0){
# log2 transform
curr_cl_data = log2(curr_cl_data)
cl_pca = prcomp(curr_cl_data)
num_cols = min(3,ncol(curr_cl_data))
curr_cl_data_raw = curr_cl_data
curr_cl_data = as.matrix(cl_pca$x[,1:num_cols],ncol=num_cols)
colnames(curr_cl_data) = paste0("clPC",1:num_cols)
for(pc in colnames(curr_cl_data)){
for(bid in curr_shared_bids){
tissue_df[tissue_df$bid==bid,pc] = curr_cl_data[bid,pc]
}
}
# Add selected CLs
for(cl in selected_cls){
if(!cl %in% colnames(curr_cl_data_raw)){next}
for(bid in curr_shared_bids){
tissue_df[tissue_df$bid==bid,cl] = curr_cl_data_raw[bid,cl]
}
}
}
}
# add sex and group, save the df
tissue_df$group = pheno[tissue_df$viallabel,"group"]
tissue_df$sex = pheno[tissue_df$viallabel,"sex"]
mt_biomarkers[[gsub("-","",tissue_name)]] = tissue_df
}
## [1] "Analyzing data from: HYPOTH"
## [1] "Analyzing data from: CORTEX"
## [1] "Analyzing data from: HIPPOC"
## [1] "Analyzing data from: ADRNL"
## [1] "Analyzing data from: LUNG"
## [1] "Analyzing data from: HEART"
## [1] "Analyzing data from: LIVER"
## [1] "Analyzing data from: SKM-GN"
## [1] "Analyzing data from: WAT-SC"
## [1] "Analyzing data from: PLASMA"
## [1] "Analyzing data from: SMLINT"
## [1] "Analyzing data from: OVARY"
## [1] "Analyzing data from: COLON"
## [1] "Analyzing data from: SPLEEN"
## [1] "Analyzing data from: SKM-VL"
## [1] "Analyzing data from: BLOOD"
## [1] "Analyzing data from: BAT"
## [1] "Analyzing data from: KIDNEY"
## [1] "Analyzing data from: TESTES"
#' Helper function: compute corr matrix from data with NAs
get_corrmat<-function(x){
n = ncol(x)
m = matrix(NA,ncol=n,nrow=n,dimnames=list(colnames(x),colnames(x)))
for(i in 1:n){
inds1 = !is.na(x[,i])
if(sum(inds1)<3){next}
for(j in 1:i){
inds2 = !is.na(x[,j])
if(sum(inds2)<3){next}
inds = inds1 & inds2
m[i,j] = cor(as.numeric(x[inds,i]),as.numeric(x[inds,j]))
m[j,i] = m[i,j]
}
}
return(m)
}
#' Helper function for the statistics
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:biomaRt':
##
## select
## The following object is masked from 'package:AnnotationDbi':
##
## select
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
compute_biomarker_dea_res<-function(x,group){
kw_p = kruskal.test(x,group)$p.value
df = data.frame(x=x,group=group)
m = lm(x~0+group,data=df)
coeffs = coefficients(m)
curr_groups = gsub("group","",names(coeffs))
curr_groups = sort(curr_groups[!grepl("control",curr_groups)])
cons = matrix(0,nrow=length(curr_groups),ncol=length(coeffs),
dimnames = list(curr_groups,names(coeffs)))
cons[,which(grepl("control",colnames(cons)))] = -1
for(r in rownames(cons)){cons[r,grepl(r,colnames(cons))]=1}
cons_res = glht(m, linfct = cons)
cons_res = summary(cons_res)
cons_coeffs = cons_res$test$coefficients
cons_tstats = cons_res$test$tstat
cons_names = names(cons_coeffs)
cons_ps = as.numeric(cons_res$test$pvalues)
names(cons_coeffs) = paste0("fc_",cons_names)
names(cons_ps) = paste0("pval_",cons_names)
names(cons_tstats) = paste0("tstat_",cons_names)
v = rep(NA,13)
names(v) = c("kw_p",paste0("fc_",c("1w","2w","4w","8w")),
paste0("pval_",c("1w","2w","4w","8w")),
paste0("tstat_",c("1w","2w","4w","8w")))
v["kw_p"] = kw_p
v[names(cons_coeffs)] = cons_coeffs
v[names(cons_ps)] = cons_ps
v[names(cons_tstats)] = cons_tstats
return(data.frame(t(v)))
}
### Stat computations
mt_biomarker_correlations = c()
mt_biomarker_dea = list()
for(tissue_name in names(mt_biomarkers)){
tissue_df = mt_biomarkers[[tissue_name]]
for(s in c("female","male")){
x = tissue_df[tissue_df$sex==s,]
biomarkers = setdiff(names(tissue_df),c("viallabel","bid","group","sex"))
for(biomarker in biomarkers){
currx = x[,biomarker]
currg = x[,"group"]
if(all(is.na(currx))){next}
inds=!is.na(currx)
currx = currx[inds];currg=currg[inds]
dea_res = compute_biomarker_dea_res(currx,currg)
dea_res$tissue = tissue_name
dea_res$sex = s
dea_res$biomarker = biomarker
mt_biomarker_dea = rbind(mt_biomarker_dea,dea_res)
}
}
curr_corrs = get_corrmat(tissue_df[,biomarkers])
rhos = c()
rhos["reads-mtdna"] = curr_corrs["mtdna","mt_reads"]
for(cl in c("clPC1","clPC2","clPC3")){
rhos[paste0("reads-",cl)] = curr_corrs["mt_reads",cl]
rhos[paste0("mtdna-",cl)] = curr_corrs["mtdna",cl]
}
rhos_cls = c()
cls = biomarkers[grepl("CL",biomarkers)]
for(clpc in c("clPC1","clPC2","clPC3")){
for(cl in cls){
rhos_cls[paste(clpc,cl,sep="-")] = curr_corrs[clpc,cl]
}
}
mt_biomarker_correlations = rbind(mt_biomarker_correlations,c(rhos,rhos_cls))
rownames(mt_biomarker_correlations)[
nrow(mt_biomarker_correlations)] = tissue_name
}
save(mt_biomarkers,
mt_biomarker_correlations,
mt_biomarker_dea,
file=paste0(data_dir,"mt_biomarkers.RData"))
Load the data instead of recomputing:
load(paste0(data_dir,"mt_biomarkers.RData"))
biomarkers = setdiff(names(mt_biomarkers[[1]]),c("viallabel","bid","group","sex"))
Fold changes and p-values:
mt_biomarker_dea[mt_biomarker_dea$tissue=="ADRNL" & mt_biomarker_dea$biomarker=="mt_reads",]
## kw_p fc_1w fc_2w fc_4w fc_8w pval_1w pval_2w pval_4w
## 13 0.24568191 -4.7687 -4.4532 -4.9286 -6.4990 0.184569199 0.1909520 0.1313492
## 15 0.01012549 9.6878 4.2518 2.4568 2.7152 0.002400526 0.2622418 0.7021477
## pval_8w tstat_1w tstat_2w tstat_4w tstat_8w tissue sex biomarker
## 13 0.03317501 -1.98126 -1.962411 -2.171908 -2.863943 ADRNL female mt_reads
## 15 0.62979979 4.01898 1.763857 1.019202 1.126400 ADRNL male mt_reads
Simple boxplots
df = mt_biomarkers$ADRNL
df$group = factor(df$group,levels = c("control","1w","2w","4w","8w"))
par(mar=c(5,4,1,1),mfrow=c(1,2))
boxplot(mt_reads~group,data=df[df$sex=="male",],las=2,xlab="",
main="Mito reads, males",ylab="Percent reads")
boxplot(mt_reads~group,data=df[df$sex=="female",],las=2,xlab="",
main="Mito reads, females",ylab="Percent reads")
for(tissue in names(mt_biomarkers)){
df = mt_biomarkers[[tissue]]
df$group = factor(df$group,levels = c("control","1w","2w","4w","8w"))
if(tissue %in% c("OVARY","TESTES","VENACAVA")){next}
num_plots = 1 # for reads
if(!all(is.na(df$mtdna))){num_plots=num_plots+1}
df$CL_72_8 = df$`CL(72:8)>CL(18:2/18:2/18:2/18:2)`
if(all(is.na(df$CL_72_8))){df$CL_72_8=df$`CL(72:8)`}
if(!all(is.na(df$CL_72_8))){num_plots=num_plots+1}
curr_col = TISSUE_COLORS[tissue]
if(tissue == "SKMGN"){curr_col = TISSUE_COLORS["SKM-GN"]}
if(tissue == "SKMVL"){curr_col = TISSUE_COLORS["SKM-VL"]}
if(tissue == "WATSC"){curr_col = TISSUE_COLORS["WAT-SC"]}
par(mar=c(5,4,1,1),mfrow=c(num_plots,2))
boxplot(mt_reads~group,data=df[df$sex=="male",],las=2,xlab="",
main=paste(tissue,"mito reads, males",sep=","),ylab="Percent reads",col=curr_col)
boxplot(mt_reads~group,data=df[df$sex=="female",],las=2,xlab="",
main=paste(tissue,"mito reads, females",sep=","),ylab="Percent reads",col=curr_col)
if(!all(is.na(df$mtdna))){
boxplot(mtdna~group,data=df[df$sex=="male",],las=2,xlab="",
main=paste(tissue,"mtDNA, males",sep=","),ylab="mtDNA level",col=curr_col)
boxplot(mtdna~group,data=df[df$sex=="female",],las=2,xlab="",
main=paste(tissue,"mtDNA, females",sep=","),ylab="mtDNA level",col=curr_col)
}
if(!all(is.na(df$CL_72_8))){
boxplot(CL_72_8~group,data=df[df$sex=="male",],las=2,xlab="",
main=paste(tissue,"CL, males",sep=","),ylab="log2 intensity",col=curr_col)
boxplot(CL_72_8~group,data=df[df$sex=="female",],las=2,xlab="",
main=paste(tissue,"CL, females",sep=","),ylab="log2 intensity",col=curr_col)
}
}
get_stat_matrix<-function(mt_biomarker_dea,stat_name,biomarkers,tissues){
m = matrix(NA,nrow=length(tissues),ncol=2*length(biomarkers),
dimnames = list(tissues,c(paste0(biomarkers,"-male"),paste0(biomarkers,"-female"))))
for(tissue in tissues){
for(biomarker in biomarkers){
for(sex in c("male","female")){
curr_ind = mt_biomarker_dea$tissue==tissue & mt_biomarker_dea$sex==sex &
mt_biomarker_dea$biomarker==biomarker
if(sum(curr_ind)!=1){next}
m[tissue,paste(biomarker,sex,sep="-")] = mt_biomarker_dea[curr_ind,stat_name]
}
}
}
return(m)
}
tissues = sort(unique(mt_biomarker_dea$tissue))
cls = biomarkers[grepl("CL",biomarkers)]
non_cls = setdiff(biomarkers,cls)
tissues = setdiff(tissues,c("OVARY","TESTES"))
ggcorrplot(mt_biomarker_correlations[tissues,c(1:7)],method = "circle")
ggcorrplot(abs(mt_biomarker_correlations[tissues,c(1:7)]),method = "circle")
ggcorrplot(mt_biomarker_correlations[tissues,-c(1:7)],method = "circle")
mt_group_kw = get_stat_matrix(mt_biomarker_dea,"kw_p",non_cls,tissues)
corrplot(t(-log10(mt_group_kw)),is.corr = F,method = "circle",
p.mat = t(mt_group_kw),sig.level = 0.05)
par(mar=c(8,8,8,8))
tissue_counts_kw = colSums(mt_group_kw<0.05,na.rm=T)
cl_inds1 = grepl("clPC",names(tissue_counts_kw)) &
grepl("female",names(tissue_counts_kw))
cl_inds2 = grepl("clPC",names(tissue_counts_kw)) & !cl_inds1
tissue_counts_kw["clPCs-male"] = sum(tissue_counts_kw[cl_inds2])
tissue_counts_kw["clPCs-female"] = sum(tissue_counts_kw[cl_inds1])
tissue_counts_kw = tissue_counts_kw[-which(cl_inds1|cl_inds2)]
barplot(sort(tissue_counts_kw,decreasing = T),las=2,col="white",
ylab = "Num tissues with p<0.05")
fc8w = get_stat_matrix(mt_biomarker_dea,"fc_8w",non_cls,tissues)
p8w = get_stat_matrix(mt_biomarker_dea,"pval_8w",non_cls,tissues)
tstats = get_stat_matrix(mt_biomarker_dea,"tstat_8w",non_cls,tissues)
corrplot(t(tstats),method = "circle",is.corr = F,p.mat = t(p8w),sig.level = 0.05)
par(mar=c(8,8,8,8))
tissue_counts_8w = colSums(p8w<0.05,na.rm=T)
cl_inds1 = grepl("clPC",names(tissue_counts_8w)) &
grepl("female",names(tissue_counts_8w))
cl_inds2 = grepl("clPC",names(tissue_counts_8w)) & !cl_inds1
tissue_counts_8w["clPCs-male"] = sum(tissue_counts_8w[cl_inds2])
tissue_counts_8w["clPCs-female"] = sum(tissue_counts_8w[cl_inds1])
tissue_counts_8w = tissue_counts_8w[-which(cl_inds1|cl_inds2)]
barplot(sort(tissue_counts_8w,decreasing = T),las=2,col="white",
ylab = "Num tissues with p<0.05")
fc1w = get_stat_matrix(mt_biomarker_dea,"fc_1w",non_cls,tissues)
p1w = get_stat_matrix(mt_biomarker_dea,"pval_1w",non_cls,tissues)
tstats1w = get_stat_matrix(mt_biomarker_dea,"tstat_1w",non_cls,tissues)
corrplot(t(tstats1w),method = "circle",is.corr = F,p.mat = t(p1w),sig.level = 0.05)
par(mar=c(8,8,8,8))
tissue_counts_1w = colSums(p1w<0.05,na.rm=T)
cl_inds1 = grepl("clPC",names(tissue_counts_1w)) &
grepl("female",names(tissue_counts_1w))
cl_inds2 = grepl("clPC",names(tissue_counts_1w)) & !cl_inds1
tissue_counts_1w["clPCs-male"] = sum(tissue_counts_1w[cl_inds2])
tissue_counts_1w["clPCs-female"] = sum(tissue_counts_1w[cl_inds1])
tissue_counts_1w = tissue_counts_1w[-which(cl_inds1|cl_inds2)]
barplot(sort(tissue_counts_1w,decreasing = T),las=2,col="white",
ylab = "Num tissues with p<0.05")
# mt_biomarker_8wfc = as.data.frame(mt_biomarker_8wfc)
# mt_biomarker_8wfc$tissue = rownames(mt_biomarker_8wfc)
# ggplot(mt_biomarker_8wfc,aes(
# x = mt_reads_female, y = mtdna_female, label=tissue
# )) + geom_point() + geom_text(hjust=-0.1, vjust=0,size=4) +
# xlim(c(-2.5,5)) + geom_abline(slope=1) + ggtitle("Female fold changes")
# ggplot(mt_biomarker_8wfc,aes(
# x = mt_reads_male, y = mtdna_male, label=tissue
# )) + geom_point() + geom_text(hjust=-0.1, vjust=0,size=4) +
# xlim(-7,6) + geom_abline(slope=1) + ggtitle("Male fold changes")
Add stats of the specific CLs
cl_tissues = c("HEART","KIDNEY","LIVER","LUNG","SKMGN","WATSC")
mt_group_kw_cls = get_stat_matrix(mt_biomarker_dea,"kw_p",cls,cl_tissues)
corrplot(t(-log10(mt_group_kw_cls)),is.corr = F,method = "circle",
p.mat = t(mt_group_kw_cls),sig.level = 0.05,main="KW test")
fc8w_cls = get_stat_matrix(mt_biomarker_dea,"fc_8w",cls,cl_tissues)
p8w_cls = get_stat_matrix(mt_biomarker_dea,"pval_8w",cls,cl_tissues)
tstats_cls = get_stat_matrix(mt_biomarker_dea,"tstat_8w",cls,cl_tissues)
corrplot(t(fc8w_cls),method = "circle",is.corr = F,p.mat = t(p8w_cls),sig.level = 0.05,
main="Week 8")
fc1w_cls = get_stat_matrix(mt_biomarker_dea,"fc_1w",cls,cl_tissues)
p1w_cls = get_stat_matrix(mt_biomarker_dea,"pval_1w",cls,cl_tissues)
tstats_cls = get_stat_matrix(mt_biomarker_dea,"tstat_1w",cls,cl_tissues)
corrplot(t(fc1w_cls),method = "circle",is.corr = F,p.mat = t(p1w_cls),sig.level = 0.05,
main="Week 1")